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ABSTRACT 

Aims. We present here a new theoretical approach to population synthesis. The aim is to predict colour magnitude 
diagrams (CMDs) for huge numbers of stars. With this method we generate synthetic CMDs for N-body simulations of 
galaxies. Sophisticated hydrodynamic N-body models of galaxies require equal quality simulations of the photometric 
properties of their stellar content. The only prerequisite for the method to work is very little information on the star 
formation and chemical enrichment histories, i.e. the age and metallicity of all star-particles as a function of time. The 
method takes into account the gap between the mass of real stars and that of the star-particles in N-body simulations, 
which best correspond to the mass of star clusters with different age and metallicity, i.e. a manifold of single stellar 
sopulations (SSP). 

Methods. The theory extends the concept of SSP to include the phase-space (position and velocity) of each star. 
Furthermore, it accelerates the building up of simulated CMD by using a database of theoretical SSPs that extends to 
all ages and metallicities of interest. Finally, it uses the concept of distribution functions to build up the CMD. The 
technique is independent of the mass resolution and the way the N-body simulation has been calculated. This allows 
us to generate CMDs for simulated stellar systems of any kind: from open clusters to globular clusters, dwarf galaxies, 
or spiral and elliptical galaxies. 

Results. The new theory is applied to an N-body simulation of a disc galaxy to test its performance and highlight its 
flexibility. 



1. Introduction 

In the past few decades the unprecedented development of 
the computational facilities drastically increased the num- 
ber of NB-based astrophysical simulations. In some areas 
of astrophysics, this powerful approach represents the only 
viable laboratory experiment because gravitational interac- 
tion on an astrophysical scale is clearly impossible to repro- 
duce in any laboratory. 

Since the pioneering works of von Hoerner (1960) and 
Aarseth (1963), the orbit integration techniques for the NB 
problem (with N the number of gravitationally interacting 
bodies) greatly improved (see, e.g., Hockney & Eastwood 
1988; Aarseth 2003) both theoretically, with the devel- 
opments of the tree-codes, particle- mesh (PM), particle- 
particle-particle mesh (P 3 M) or see, e.g., Dehnen & Read 
(2011) for a recent review, and technically, by exploiting the 
message-passing-interface protocol for the use of massively 
parallel machines. Although that the N-body simulations 
were originally conceived as a tool to follow the integration 
of the orbits, the importance of including energy dissipative 
processes soon became evident. The treatment of baryonic 
interactions in N-body simulations is nowadays standard- 
ised by popular protocols such as smoothed particle hy- 
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drodynamics (SPH) or adaptive-mesh-refinement 1 (AMR) 
implemented in several codes, e.g., Algodoo (e.g., Kores 
2012), DualSPHysics (e.g., Gomez- Gesteira et al. 2012), 
SPH-ftow (e.g., Oger et al. 2006), ENZO (e.g., O'Shea 
ct al. 2004), EvoL (e.g., Merlin et al. 2010), GCD+ (e.g., 
Kawata k Gibson 2003), GADGET (e.g., Springel et al. 
2001), Gasoline (e.g., Wadsley et al. 2004), FLASH ? or 
RAMSES (e.g., Teyssier 2010; Few et al. 2012), which 
also includes (in astrophysical context) a sticky-particle ap- 
proach (Bournaud & Combes 2002; Martig & Bournaud 
2010). Thus, even though the recipes for implementing the 
physics of dissipative phenomena are still controversial, sev- 
eral NB codes are currently able to evolve with time the 
baryonic component, i.e. gas and stars, in their mutual in- 
teraction. Therefore, it is nowadays possible and manda- 
tory to develop the right tools to compare the results of 
dynamical N-body simulations with observational data for 
the stellar content of galaxies. 

Tantalo et al. (2010) recently developed a technique to 
derive the integrated spectra, magnitudes and colours of 
the stellar content of simulated galaxies. These authors fo- 
cused on the evolution of non-resolved stellar populations 
and exploited the concept of the spectral energy distribu- 



1 Although AMR codes are not an NB-type code, they often 
use star-particles to describe the stellar components. Therefore, 
we here include AMR for our terminology of N-body simulations 
in this paper. 
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tion (SED) of an SSP in the context of discrete resolved- 
mass points of an N-body simulation. There, the integrated 
monochromatic flux in a given photometric system (set of 
discrete pass-bands AA) of a galactic stellar component 
at the time t and metallicity Z, F\ (Z;t), is taken to be 
the convolution of the star formation rate (SFR) tp(Z;t) 
(which is assumed to depend only on age and metal con- 
tent), and the integrated monochromatic flux of constituent 
SSPs, 4> x (Z; t), i.e. F x (Z; t) = V (Z; t) * 4> x (Z; t) 2 . 

In this study we present a new theory of population 
synthesis that is particularly designed to manage very 
many stars, and is based on the concept of a distribu- 
tion function (DF) . We focus then on generating the CMDs 
for the stellar system which are simulated with advanced 
hydro-dynamical N-body simulation techniques. The sub- 
ject is particularly timely in view of the modern capabilities 
of data acquisition of star-by-star photometry in nearby 
galaxies. This is thanks to modern instrumentation already 
in place or in progress for the coming years, and also in 
view of the N-body simulations of model galaxies, in which 
sampling of either real or fake stars with many millions of 
objects is possible. 

The only prerequisite for the method to work is that 
minimal information about the history of star formation 
and chemical enrichment is implemented in and read off 
the N-body simulations. For each star-particle of the N- 
body simulation we must know the age at which it was born 
and the chemical composition (metallicity) of the parental 
medium. The new formulation of population synthesis takes 
naturally into account the gap between the mass of the real 
stars and the mass of the star-particles in N-body simula- 
tions if these latter are closer to the mass of a star cluster 
than the mass of a real star. Therefore the star-content of 
an N-body simulation is better described as a manifold of 
star clusters with different age and metallicity, i.e. a man- 
ifold of SSP whose photometric properties, primarily the 
CMD, are well known. Second, CMDs containing an arbi- 
trary number of stars can be easily simulated by defining 
and making use of the concept of DF of stars in the CMD. 

This novel technique is applicable to the stellar con- 
tent of single open/globular clusters, to our own Galaxy, to 
that of nearby galaxies within the Local Group, and to all 
galaxies of the local Universe whose stellar content can be 
resolved into stars. Moreover, the algorithm can be easily 
interfaced with other codes, e.g. Galaxia (e.g., Sharma et al. 
2011), or the Galaxy Simulators HRD-GST (e.g., Ng et al. 
1995, 1996; Ng & Bertelli 1996; Ng et al. 2002a,b; Vallenari 
et al. 2006), and Besancon model (e.g., Robin et al. 2003), 
thus extending their performances with the advantages that 
will become soon clear in what follows. 

Setting up the new technique, we focused our attention 
on a few key requirements that should be met: 

1. The algorithm has to be independent of the particular 
N-body simulation it is applied to. Thus, it is designed 
to accept the output of a simulation as input, without 
interfering with the NB integration. 



2 Where / * g is the standard convolution operator between 
the generic functions / and g. In general, we can indeed always 
write F\ = L^ 1 [L [tp] L [4>\]] , with L [ip] the discrete Laplace 
transform of the SFR we obtain directly from the N-body simu- 
lations and L [<j>\] the Laplace transform of the monochromatic 
integrated flux of the a simple stellar population. 



2. The algorithm must not depend on the specific pre- 
scriptions used to generate the stellar models that are 
adopted to describe the composite stellar populations 
and associated CMDs. 

3. The technique needs to handle CMDs populated by sev- 
eral billions of stars, i.e. which are potentially represen- 
tative of the largest systems of stars known (e.g., giant 
elliptical galaxies) . Indeed, our algorithm should be ap- 
plicable to every resolved stellar population no matter 
what the nature of the stellar system under considera- 
tion. Therefore CMDs for giant elliptical or spiral galax- 
ies that contain up to 10 12 stars have to be synthesized 
with agility. 

4. Finally, we require the code to only require little compu- 
tational resources. The typical CMD of point 3, has to 
be obtained on a serial machine, i.e. its realization has 
to avoid more complicated MPI/OpcnMP programming 
protocols. 

To achieve those goals, we have developed a method for 
synthesising CMDs that takes all the above requirements 
into account. The plan of the paper is as follows: in Section 
2 we extend some of the theory concepts of the stellar pop- 
ulations within the framework of the DFs. In Section 3 we 
particularise the general theory to the case of the N-body 
simulations. In Section 4 we deal with the simulations of 
CMDs that contains huge numbers of stars. In Section 5 
we present an example of a synthetic CMD realization. In 
Section 6 we show an example of the technique. Finally 
we draw some concluding remarks in Section 7 and briefly 
mention some possible applications of the new method to 
the synthesis of stellar populations. 

2. Theory of Stellar Populations 

We extend here the theory of stellar populations to systems 
(assemblies of stars) with an assigned distribution in the 
phase-space. This completes the standard theory of stellar 
populations (e.g., Salaris & Cassisi 2005; Greggio & Renzini 
2011) by assigning a position and velocity to each star. 

Every real (or realistically simulated) stellar system is 
a set of stars born at different times and positions, and 
with different velocities, masses and chemical compositions. 
We call this assembly of stars composite- stellar- population 
(CSP). The position x of each star of a CSP is a point 
in the space X of coordinates (or space of configurations), 
and its linear momentum p is geometrically referred to 
as a fiber of the cotangent space T*X of X at x (Jose 
& Saletan 1998). Therefore, the phase-space of a CSP is 
the cotangent bundle T*X = 1JT*X of its configuration 

X 

space X. It follows that T*X is a dim[T*X] = 67V mani- 
fold, with N the total number of stars in the CSP. Each 
point in it defines the phase-space of our CSP. We indi- 
cate this manifold with the symbol Y, i.e. Y = T*X where 
T = (x,v) = (x 1 ,x 2 ,...,x 3N ,v 1 ,v 2 ,-.,v 3N ). Moreover, to 
completely define the CSP parameters at a generic time 
t, we need to specify the distribution in the space of the 
masses, M, and metallicities, Z of all its members. We call 
this extended phase-space the existence- space of the CSP, 
i.e. E = M x Z x Y of the CSP at the time t and Exl 
the extended- existence- space with time t. In E x R, the stars 
evolve with time, i.e. they continuously move in space, lose 
mass, enrich in metals, and move in the phase-space. Hence 
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we can safely define in E the distribution function (DF) for 
the CSPs under the assumption of continuity and differ- 
entiability. Let us consider a sample of identical systems 
whose initial conditions span a certain volume of the space 
E. We refer to this sample simply as an ensemble, borrow- 
ing the name-root from the grand microcanonical ensemble 
adopted in Statistical Mechanics. The number of systems 
dN at the time t with mass within dM, metallicity within 
dZ and phase-space within dT = (dx, dv) is given by 



dN = Nfcsp (M, Z, T; t) dMdZdT, 
with /csp (M, Z, T; t) the DF in E, / CS p 



(1) 



p6W 



. X K X 

— > R + continuous and with partial derivative 
continuous (where R is the set of all real numbers, and 
R + = {a\a 6 MAo ^ 0}) and the total number of systems 
in the ensemble is fixed by normalising to one the DF in 
the space E, i.e. 



J /csp (M, Z,T;t) dMdZdT 



1. 



(2) 



We proceed to formally define an SSP as follows. At a given 
time t, we divide E into a grid in terms of discrete intervals 
of metallicity dZ and phase-space dr. Then every infinites- 
imal unit (i.e. every elementary cell) of this sub-space of E 
defines an SSP. 

We assume that every CSP originates at time to from 
an episode of single star formation in a chemically homoge- 
neous medium. The stars of the CSP have the mass spec- 
trum f M = f M (M,t Q ) oc f (M,t ) where £(M,t ) is the 
initial mass function (IMF). Their metallicity distribution 
is fz = fz (to) cx Zqsp {to)- Finally, their phase-space dis- 
tribution function is / r = /r (r;i ). The DF of the CSP 
can be written as 



/csp = X] & 



SSP, 



(3) 



where n is the number of SSP, with eventually n — > oo, 
and / SSP = /ssp (M, Z , T ; t ) is the DF of SSP of which 
Z = Z , r = r , t = to are the metallicity, SSP's phase- 
space and age at the instant to, respectively. 

As time elapses, this stellar population ages. According 
to their mass, the stars eventually leave the main-sequence 
(MS) after a time tyis = tMs(M) and soon after die (super- 
novae) or enter a quiescence stage (white dwarfs), inject- 
ing metals into the interstellar medium in form of super- 
nova explosions or quiet winds. The interstellar medium 
becomes richer in metals due both to self-enrichment by 
the SSP under consideration and the contribution from all 
other stellar SSPs. At the generic time t > to, i.e. after an 
elapsed time r = t — to otherwise known as the age of the 
stellar population, the SSP has processed stellar material. 
As a result of this activity, an age-metallicity relationship 
Z CS p = Zcsp (*) is built up for the CSP. 

In addition to this, the number of trajectories in the 
phase-space entering a volume dT, will in general be dif- 
ferent from the number of those leaving the same volume. 
Thus /r evolves following Liouville's equation of the form 



dfr 

dt 



d_ 

or 



r-^)/r, 



(4) 



which is independent of the nature of the equation of mo- 
tion (e.g., its correctness in the form of Eq. (4) does not re- 
quire the existence of a Hamiltonian). C is the Liouvillean 



operator, eventually extended to account for a creation 
function C (/csp;*), an d i the imaginary unit. The formal 
solution of this equation is given by the Taylor expansion 
series of the time dependence of /r (T; t) around /r (T; t ) 

fr = e-^/r (T;t ) = £ ^j-^/r (T;t ), (5) 

m=0 



and 



-iCt _ 



E 

m— 



(-*)" 



■(tC) m defines the infinite series of 



operators applied on any function on its right side. In this 
way, the initial SSP has generated a CSP whose distri- 
bution in the existence space of the CSP, E, has a DF 
f CS p(M,Z,T;t). 

Now that we have stated these fundamental concepts, it 
is easy to proceed with the following formal definitions to 
recover the usual concepts of the classical synthesis of pop- 
ulation theory: 

— Present-day-mass-function: The integral of the DF over 
the metallicity Z and phase-space T (from Eq. (I)) 



/ 



iV/csp (M, Z, T, t) dZdT = £ (M; t) , 



(6) 



yields the "present" -day- mass- function £(M;t) 
(PDMF), i.e. the total number of stars per mass 
interval at the time t. This can be expressed by the 
approximate relation 



i(M;t) = 



£(M) 



*MS 
t-t 



tMS < T 



£ (M) t MS > T, 



(7) 



where £ (M) is IMF of the MS stars. More fine-tuned 
approaches based on the evolutionary flux, total 
luminosity and fuel-consumption theorem can be easily 
worked out upon necessity. 

Age-metallicity function: Integrating the DF over the 
mass M and phase-space T 



I 



Nfcsp{M,Z,T;t)dTdM = X {Z,t), (8) 



we obtain the age-metallicity relation, i.e. the num- 
ber of stars formed per metallicity interval at the time t. 

The phase-space distribution-function: Integrating upon 
the mass M and metallicity Z 



J M/csp (M, Z, T, t) dMdZ = e- LCt f T (T; t ) , 



(9) 



where we made use of Eq. (5), yields the the phase-space 
DF. 

In a natural way, we can extend the approach used in 
defining Eqs. (6), (8) and (9) by performing the following 
monodimcnsional integration: 

Metallicity phase-space relationship: We consider the fol- 
lowing integration 



(10) 



V (Z,T;t) = / f CS p(M,Z,T;t)dM 
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to define the metallicity-phase-space relationship. Very of- 
ten the projection of this function onto the configuration 
space is used 

J V (Z,T;t)d m v = f,(Z, X ;t) 7 (11) 

where d 3N v is the 3N velocity element of the T space. This 
plays an important role, e.g., in relation to radial metal- 
licity gradients of resolved populations in external galaxies 
and the vertical metallicity gradients in the Milky Way, i.e. 
fj (Z, x;t) — Vf can implicitly define a function that re- 
lates the metallicity Z and the radius of a galaxy R once 
cylindrical coordinates are in use. 

Mass-phase-space relationship: In the same way we define 

H (M, T;t)= J / CSP (M, Z, T; t) dZ, (12) 

i.e. the mass-phase-space relation. More often its projection 
onto the configuration space 

J ji (Z, T; t) d m v = pi (Z, x; t) (13) 

is used, which is related, e.g., to the evolution of the globu- 
lar clusters and the mass-segregation effects, i.e. we can use 
it to express the higher concentration of massive (or binary 
stars) through the centre of the cluster with respect to the 
less massive stars (e.g., Spitzer 1987). 

Mass-Metallicity relationship: Finally, we define 

w(M, Z;t)= J ic S p (M, Z, T, t) dT, (14) 

i.e. the mass-metallicity relation for a resolved stellar popu- 
lation. This relation is of paramount importance e.g., when 
the masses of the stars are obtained by parallaxes, transits, 
and astero-seismology. For nearby galaxies this relation is 
relevant only in studies of integrated photometry. 

3. Stellar populations in N-body simulations 

The concepts we have just developed can be straightfor- 
wardly applied to N-body simulations of any kind, from 
those for globular clusters with single or multiple stellar 
populations (e.g., Milone et al. 2012; Norris 2004), to those 
for the Milky Way stellar fields observed along any line 
of sight (e.g., Ng et al. 1995; Vallcnari ct al. 2006), to 
the dwarf galaxies of the Local Group (e.g., Grebel 1997; 
Mateo 1998; Tolstoy et al. 2009), or even external galaxies 
once they are resolved into individual stars (e.g., Harris & 
Harris 2002; Rejkuba et al. 2011; Crnojevic et al. 2011). 
Indeed, in the adopted formulation, all the dynamical 
aspects of the processes governing the formation of the 
astronomical object under investigation are considered 
through the most general Liouvillean operator, whose form 
can be specified case by case. In this context, it is worth 
recalling that a mutual dependence of star formation, 
which is responsible for the building up of the stellar pop- 
ulation of a stellar system, and the dynamical processes 
that govern the large scale aggregation of it, is in general 
supposed to occur in extended many-body interactions 
(e.g. Pasetto et al. 2011, 2012). 



Star-particles in N-body simulations: Before proceeding 
further, we must suitably link the definition and proper- 
ties of CSPs to the elemental building blocks of N-body 
simulations, i.e. the "star-particles". 



We define an N-body simulation MS M stochastic real- 
ization of /csp (M, Z,T;t) in the existence space E of the 
CSPs. We focus our attention on the evolution of a sin- 
gle star-particle in an N-body simulation. This particle p is 
born at the instant to, p = t p with metallicity Z p = Z p (t p ), 
at a position in the single particle phase-space 7 (a sub- 
manifold of T*X) given by 7 P = 7 P (i p ), with a certain 
metallicity Z = Z(t p ) (the metal content of the parent gas 
component at time t p ), and with a mass M p = M p (t p ). 



The mass of the star-particles and the smaller mass res- 
olution M of the N-body simulations is one of the great- 
est problems of modern computational N-body simulations 
and companion mathematical techniques. Apart from a 
few recent N-body simulations that were specifically tai- 
lored for star clusters (e.g., Hut et al. 2003; Zonoozi et al. 
2011; Pelupessy & Portegies Zwart 2012), the particle mass 
of the modern N-body simulations dedicated galaxies and 
galaxy clusters can vary within a wide range, typically 

log 10 (m) -2^-8 (e.g., Saitoh et al. 2009, 2008; Guedes 

et al. 2011; Okamoto & Frcnk 2009). This means that the 
average mass of the star-particles in the NB simulations is 
larger than the average mass of the stars in a SSP. Therefore 
each star-particle is representative (in mass) of many hun- 
dreds stars (or more). However, since all stars in a star- 
particle of a simulation are assumed to be born at the same 
time with the same metallicity, the stellar content of the 
star-particle itself is well described by an SSP (see SSP def- 
inition in Section 2). Nevertheless, NB simulations cannot 
easily describe the phase-space distribution of individual 
stars within a star-particle. In this paper we assume that 
all stars within each star-particle share the same position in 
the phase-space. As a consequence of this, the volume of the 
existence space E occupied by elemental cells is bounded by 
the mass resolution, say dE m ; n = dM p x dZ x dF, and it 
is generally larger than the volume necessary to properly 
map a CSP, dE = dM x dZ x dT. 



On one hand, if the problem clearly shows that attempts 
to map the complete IMF are at present out of reach, this 
indicates on the other hand how a possible solution has to 
be searched for in the concept of SSP itself. 



In an N-body simulation, at the time t > t p a star- 
particle is located in -f p — j p (t), has metallicity Z p = Z p (t) 
and mass M p — M p (t). Therefore, mass, metallicity and 
position in the single-particle phase-space univocally iden- 
tify the NB-particle all along its evolutionary history from 
t p = ta. P to t. At the initial time t p for each particle p we 
can associate an SSP to the star-particle in a natural way, 
provided that some re-scaling of the SSP mass, Mssp, is 
made in relation to the star-particle mass, M p . From Eq. 
(2), the total mass of the SSP can be written as follows: 
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MNf CSP (M, Z,T;t)8(Z 
= J MNJcsp (M,Z p ,T p ;t p )dM 
= J Mi(M;t p )dM 
= J MZ(M)dM = M SSP (t p ), 



z p , r 



T p ,t-t p )dMdZdT 



where S is the ordinary multidimensional Dirac function, 
and in the last row we used Eq. (7) with t p = io,p — < 
tMsVM. In other words, the SSP mass written in terms of 
the DF /csp or the IMF. 

In the same way, we can cast the DF /csp (M, Z, T; t) 
of the CSP as a function of the star-particle mass M p and 
vice versa. The following identities can be written: 



/ 



M r , 



-MNfcsp (M, Z p , T p ] t p ) dM 

M S SP 

= J MNf CS p, P (M t p )dM 

= j MNf SS p } p(M,Zp,r p ;tp)dM 
= M p (t p ) . 



(16) 



While associating the SSP to a generic star-particle born 
at time t p with metallicity Z p (t p ) = Z p and mass M p (t p ) — 
M p is straightforward, at a generic time t > t p it is no 
longer so simple. But if the mass and metallicity of a star- 
particle do not change with time, the SSP associated to a 
generic star-particle will simply evolve and age, losing a cer- 
tain amount of mass in form of gas, changing a fraction of 
its living stars into remnants, and changing its spectral and 
photometric properties in a way that is fairly well known 
from the theory of stellar evolution and population synthe- 
sis (e.g., Tosi et al. 1991; Bertelli et al. 1994; Tolstoy & 
Saha 1996; Aparicio et al. 1997; Harris & Zaritsky 2001; 
Gallart et al. 2005; Bertelli et al. 2008, 2009). 



4. Huge CMDs at no cost 

With the formalism we developed on the basis of the CMD 
of CSP, we now tackle the key question: how can we build 
CMDs for billions of stars in practice at less computational 
cost? Because real CPS are the result of the past history 
of star formation in clusters and fields at the complexity 
levels of galaxies, how can we easily achieve the same level 
of information in numerical simulations? 



4.1. Associating an SSP to a star-particle 

The SSP corresponding to the p th star-particle is obtained 
by two-dimensional interpolation (in age and metallicity). 
The number of stars dv v (L, T c g) for the interval of luminos- 
ity L and effective temperature T e g, dLdT e g results from 



where we referred, by extension of the Jacobean-matrix for- 
malism, with 



d(r,Z) 



_d(T e{{ ,L) 

two of the dimensions in which SSP is defined 



to the transformation applied on 
from age- 

metallicity space to the effective temperature - luminosity 
space. Or analogously, 



dvp (C,m) = / SSP , p [[j^r^W [[-d(^ny\\ dCdm 

(18) 

for the same transformation to the space of the colour C 
and magnitude m. In general there is no analytic formula- 

They 



tion for the two matrixes 



d(r,Z) 



and 



d(T e it,L) 



d(C,m) 



_[d(T e!! ,L)_ 

are derived numerically from the tabulations of bolometric 
corrections, and colours of, e.g., the Johnson-Cousin-Glass 
system. 

The key idea of Eq. (17) or (18) is to describe a CMD as 
a matrix, i.e. a projected DF, whose elements are the rel- 
ative frequency (or percentage) of stars of different colour 
and magnitude in some photometric system, per elemental 
area of the CMD. For a given photometric system with pass- 
bands AA Q (for instance a stands for U, B, V,...) for which 
magnitudes m a and colours C Q( g = m a — mp can be com- 
puted (magnitudes and colours can be either absolute and 
intrinsic or apparent and reddened 3 ), Eq. (18) represents a 



2D grid of elemental square cells dv p (C°p,m^j , identified 

by the coordinates of their centres m c a and (m a —mp) c . The 
path drawn in the CMD by a single SSP of assigned age 
and chemical composition, see Section 2 and Eq. (18), will 
occupy a number of cells from the main sequence to the last 
observable stage. Each cell is populated by a certain num- 
ber of stars according to the underlying luminosity function 
of the SSP, which in turn is related to the evolutionary rate 
and IMF (see Eq. (6) and (7)). If many SSPs are present, 
as in the case of a composite CSP, each cell contains a to- 
tal number of stars given by the sum of all contributions 
by different SSPs passing through this cell. Using Eqs. (3) 
together with (18), we obtain 



Vp (C a0 ,m a ) = ^dv p (C a p,m a ) 



(19) 



Therefore, a relative percentage of stars is associated to 
each cell (matrix element) with respect to the total. The 
number of stars per cell can be easily normalised accord- 
ing to the problem under investigation. The observational 
CMDs that contain an arbitrary number of stars, from a 
few thousands to billions, is reduced to a matrix of number 
frequencies (relative percentages, Eq. (18)). The prescrip- 
tion can be easily extended to include any history of star 
formation of any intensity, ip, as well as kinematical effects 
(fr as in Eq. (5)). 

Finally we are left to specify the initial distribution of 
mass (IMF) and the adopted database of stellar models: 

— The initial mass function: To calculate an SSP, we need 
an IMF. The most popular IMF is given by 



^(M) = ^Y, M ~ Xi ^ 



(20) 



dvp (T e ff, L) = /ssp.p 



d(r, Z) 



d(T cS ,L) 



dT cS dL, 



(17) 



3 The same considerations also apply to bolometric luminosi- 
ties and effective temperatures. 
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Fig. 1. Left panel: The HRD of SSPs with assigned age and metallicity. The size of the circles indicates the current star mass along 
the isochrones and/or SSP. The stellar models in use are taken from the library of Bertelli et al. (2008) and Bertelli et al. (2009). 
The colour code refers to the size of the dots: the smaller is the dot, the more yellow it is and the smaller the mass of the plotted 
stars. The bigger the dot, the bluer its colour and the bigger the mass of the stars. The mass considered in the models of this 
example ranges from 0.15 to 20 M Q . Right panel: The same SSPs displayed in the left panel but in the observational CMD My vs 
My — Mi. The CMD is subdivided to a discrete grid of elemental cells, in each of which a certain number of stars fall. The stars 
in each cell may belong to different SSPs. The size of the cell is large clearly show evidence the "out-of-focus" effect intrinsic to 
this tessellation technique. Colour-code spans account for the number of stars in the cells from (the bluer colour) to 10 5 (the 
light yellow colour). 



with 



1 M G AMi 
M <£ AMi 



(21) 



with £ normalisation factor of the IMF, AM, = 
{Mi\Mi G [Mi ow ,i,M up ,i\}, with M JoW) < (M uPji ) lower 
(upper) limit for the i -mass interval and AM; n 
AM j = for i ^ j. The case k = 1 in Eq. (20) is of- 
ten referred to as Salpeter's IMF with power-law slope 
xi = 2.35 Vi, the case k = 3 is referred as Kroupa's IMF 
with its specific slope and intervals (see e.g., Kroupa 
2001) or in general any IMF function can be approxi- 
mated by a sum as in Eq. (20) for a suitable choice of k. 
Therefore, in the general case, the total number of stars 
plotted in an SSP is 



SSP 



AI n 



Zo^AiM-^dM 



i=l 



M~ Xi dM 



(22) 



M,, 



k f M \ 

or N SS p = £ J2 log ( for Xi = 1. If we choose 

i— 1 ^ ° w ^ 

to sample all star-phases of higher luminosity (e.g., 



asymptotic giant branch (AGB) and planetary nebula 
(PNs)) with iVssp — 10 5 stars it becomes immediately 
evident how the graphical realization of a CMD for an 
N-body simulation encounters serious problems once 
the number of star-particles involved, N, is high, say 
iV SSP x N = O (10 11 " 13 ) for N - 10 6 ^ 8 . 

The database of SSPs: We briefly report here on the data 
base of SSPs that was calculated for the purposes of 
this study. The stellar models used are those by Bertelli 
et al. (2008, 2009), which cover a wide grid of helium 
Y, metallicity Z, and enrichment ratio AY/ AZ. The 
associated isochrones include the effect of mass loss by 
stellar wind and the thermally pulsing AGB phase ac- 
cording to the models calculated by Marigo & Girardi 
(2007). The code used is the last version of YZVAR 
developed over the years by the Padova group used 
in many studies (for instance Chiosi & Greggio 1981; 
Chiosi et al. 1986b,a, 1989; Bertelli et al. 1995; Ng 
et al. 1995; Aparicio et al. 1996; Bertelli & Nasi 2001; 
Bertelli et al. 2003) and was recently extended to ob- 
tain isochrones and SSPs in a large region of the Z — Y 
plane. The details on the interpolation scheme at a given 
AY/AZ are given in Bertelli et al. (2008, 2009). 
The present isochrones and SSPs are in the Johnson- 
Cousins-Glass system as defined by Bessell (1990) and 
Bessell & Brett (1988). The formalism adopted to derive 
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the bolometric corrections is described in Girardi et al. 
(2002), while the definition and values of the zero-points 
are described in Marigo & Girardi (2007) and Girardi 
et al. (2007) and will not be repeated here. 
Suffice it to recall that the bolometric corrections stand 
on an updated and extended library of stellar spec- 
tral fluxes. The core of the library now consists of 
the ODFNEW ATLAS9 spectral fluxes from Castelli 
& Kurucz (2003), for T off G [3500,50000] K, log lo5 £ 
[—2,5] (with g the surface gravity), and scaled solar 
metallicities [M/H] e [-2.5, +0.5]. This library is ex- 
tended at the intervals of high T c g with pure black-body 
spectra. For lower T c g, the library is completed with 
the spectral fluxes for M, L and T dwarfs from Allard 
et al. (2000), M giants from Fluks et al. (1994), and fi- 
nally the C star spectra from Loidl et al. (2001). Details 
about the implementation of this library, and in partic- 
ular about the C star spectra, are provided in Marigo & 
Girardi (2007). It is also worth mentioning that in the 
isochrones we applied the bolometric corrections derived 
from this library without making any correction for the 
enhanced He content which has been proved by Girardi 
et al. (2007) to be low in most common cases. 
The database of SSP covers the space of exis- 
tence E of a generic CPS. The number of ages 
N T of the SSPs are sampled according to a law 
of the type r = i x 1CP for i = 1,...,9 and 
j = 7,..., 9, and for Nz metallicities are Z = 
{0.0001, 0.0004, 0.0040, 0.0080, 0.0200, 0.0300, 0.0400}. 
The helium content associated to each choice of 
metallicity is according to the enrichment law 
AYAZ = 2.5. Each SSP was calculated allowing 
a small age range around the current value of age 
given by At = 0.002 x 10-? with j = 7, ...9. In total, 
the data base contains N T x N z = 150 SSP. This grid 
is fully sufficient to illustrate the method. For future 
practical application of it, finer grids of SSPs can be 
calculated and made available. Having done this, the 
normalization constant £ of the SSPs remains defined. 
Finally, for each SSP we computed the "projected" DF, 
i.e. the number of stars per elemental cell of the CMD, 
using the value of £o- 

Note that the method for generating the cumulative 
DFs of Eq. (19) from SSPs does not depend on the 
particular choice for the data base of stellar tracks, 
isochrones, or photometric system. Other libraries of 
stellar models and isochrones can be used to generate 
the database of SSP, the building blocks of our method. 

From the procedure outlined above, it is easy to under- 
stand how the use of the stellar DF is able to accelerate 
the construction of the CMD of a CSP. The reason is as 
follows: Instead of calculating an SSP made by iVssp stars 
for every star particle of the N-body simulation (i.e., for 
a total of N star particles) and counting the stars inside 
a given bin of magnitude and colour 5C5m, i.e. summing 
over N x Assp — O (lO 11 ^ 13 ) stars, with the above proce- 
dure, we need to calculate only N T x Nz x -/Vgsp = O (10 7 ) 
stars. However, the number frequencies per cell of the CMD 
of our reference SSPs are calculated once for all, whereas 
their combinations can be freely changed according to the 
underlying star formation history of the N-body simulation 
to investigate and the number of calculation required in this 
method is then just O (N) < 10 11 " 13 . 



The left panel of Fig.l shows a few SSPs for the solar 
metallicity Z = 0.02 and helium content Y = 0.28 in the 
theoretical Hertzsprung- Russell diagram (HRD). The size 
of the dots is proportional to the stellar mass running along 
the SSPs. The same SSPs are translated into the My vs 
My — Mj plane displayed in the right panel of Fig. 1 in 
which the cell tessellation is evident. The "out-of- focus" 
effect is due to grouping the stars of different SSPs into the 
same cell. No photometric errors are applied. 

4.2. Simulation of photometric errors and completeness 

We finally mention that real data on the magnitudes (and 
colours) of the stars are affected by photometric errors, 
whose amplitude in general increases at decreasing lumi- 
nosities (increasing magnitude). The photometric errors 
come together with the data themselves provided they are 
suitably reduced and calibrated. Photometric errors can be 
easily simulated in our theoretical CMDs. Suppose we know 
the errors affecting our magnitudes in the two pass-bands 
used to build the CMD, and that these are a function of 
the magnitude itself. Let us indicate the errors by Sm a (m), 
and 5mp(m), with a and (3 the two pass-bands. When plot- 
ting each point of the SSPs on the observational CMDs, the 
magnitudes are changed by the quantity representing the 
errors, e.g.: 

m' a (m) = m a (m) - Q - r^j 6m a (m) 

m!p{m) = mp{m) — ^— — s^j 5mp{m), (23) 

where r and s are two randomly drawn numbers (e.g., from 
uniform or Gaussian distribution) comprised between and 
1. The star frequencies per elemental cell of the CMD are 
calculated after applying the correction for photometric er- 
rors to the reference SSPs. 

To compare the real observational data with the theoret- 
ical model we have to know the completeness of the former 
as a function of the magnitudes and pass-band (Stetson & 
Harris 1988; Aparicio & Gallart 1995). This is a long known 
problem that docs not require any particular discussion in 
the context of this paper, and tabulations of the complete- 
ness factors must be supplied in advance together with the 
correction for photometric errors. We mention here is that 
correcting for completeness will alter the DF of stars in 
the cells of the observational CMD we aim to analyse (e.g., 
Crnojevic et al. 2011). These completeness factors must be 
supplied by the user of our method in connection with the 
specified problem. 

5. Generating a CMD with tessellation 

To explain the concept of CMD tessellation, we simulated 
a CMD using the CSP of NB-TSPH simulation (to be de- 
scribed in more detail below). The analysed field contains 
8.8 x 10 4 star-particles of the same mass and known age and 
metallicity (to each of which an SSP can be associated with 
the same mass, age and metallicity, see Section 3). For the 
purpose of this example we show only the absolute luminos- 
ity and magnitude, i.e. no distance dependence. The results 
of applying Eq. (19) are shown in Fig. 2. The left panel dis- 
plays the frequency distribution of CSP in the (My — Mj) 
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Fig. 2. Left panel: The CSP for an N-body simulation of a disc galaxy. Colour code refers to the frequency of star per bin shown in 
the companion right panel. Right panel: The histogram of DF for the CSP shown in the left panel. In this figure, one immediately 
captures the concept of DF described in the text. The characteristic peaks corresponding to MS and RGB stars are shown in 
yellow. The CMDs have the coordinates of the absolute magnitudes Mv and the colours My — Mi. 



$ ° 



■ft 
I 




Fig. 3. Snapshot of the simulated disc galaxy whose face-on views of the stellar and gas discs are shown in the left and right 
panels, respectively. Our assumed location of the Sun is highlighted by a star symbol. The white lines indicate the selected region 
for constructing the CMD shown in Fig. 5. Colour code ranges from blue (the lowest density) to red (the highest density). 



vs My CMD, while the right panel shows the histogram. 
The regions occupied by stars in the main sequence are 
clearly evident: Long-lived phases display a higher num- 
ber of stars, on the other hand, red giant, red clump, and 
asymptotic giant phases are seen in low effective tempera- 
tures (red colours) and shows lower frequency due to their 
short lifetime. These features are similar to a typical ob- 
served CMD of the stellar populations in nearby galaxies, 



e.g. the Magellanic Clouds, M31 and others, where for all 
stars in a given galaxy the distance is nearly constant. 

6. Application to N-body simulations 

To demonstrate the method, we applied it to an N-body 
simulation. The simulation was carried out with an updated 
version of our original NB-SPH code, GCD+ (Kawata & 
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Gibson 2003; Rahimi & Kawata 2012). We initially set up 
an isolated Milky Way-sized disc galaxy that consists of gas 
and stellar discs with no bulge component in a static dark 
matter halo potential, following Rahimi & Kawata (2012). 
Note that this simulation was used for demonstration pur- 
pose, and was not meant to reproduce the Milky Way. 
We used the standard Navarro-Frenk- White dark matter 
halo profile (Navarro et al. 1997) with the total mass of 
Mtot = 1.5 x 10 12 Mq and the concentration parameter of 
c = 12. The mass, scale length and scale high of the stellar 
disc were assumed to be M d>s = 4.0 x 10 10 M Q , i? djS = 2.5 
kpc and z djS = 350 pc. The mass and scale length of the gas 
disc was M d>g = 1.0 x 10 10 M Q , i? d , g = 4.0 kpc. We initially 
set 400,000 particles to the gas disc and 1,600,000 particles 
to the stellar disc. Therefore, our baryon particle mass was 
M p = 2.5 x 10 4 M Q . We applied the threshold density of 
7iH = 1-0 cm -3 for star formation. Rahimi & Kawata (2012) 
demonstrated that the star formation in a disc is quite sen- 
sitive to the parameters of star formation efficiency, C*, 
energy for supernova, -Esn, and stellar wind feedback en- 
ergy, Esw- We here applied C* =0.1, = 10 51 er g an d 
Esw = 10 37 erg s _1 . Fig. 3 shows a snapshot of the sim- 
ulation. Strong feedback creates many bubbles in the gas 
disc. 

The simulated galaxy shows a small bar and several spi- 
ral arms. We set the Sun at (x, y) = (—8, 0) kpc (star sym- 
bol in Fig. 3), and selected star-particles in the region of 
the simulated galaxy "equivalent longitude" 300 < I < 320 
and "latitude" — 10 < b < 10 (with implicit reference to 
the Milky Way galaxy coordinates system I, b), which is en- 
closed by white lines in Figs. 3 and 4. For each star-particle 
we measured the distance and extinction from the simu- 
lation. To measure the extinction, first the column den- 
sity for each star-particle was calculated by summing up 
the line of sight column density of all gas particles be- 
tween the star-particle and the position of the Sun, using 
the SPH weighting scheme (Kawata & Rauch 2007). We 
then converted the column density to the extinction, using 
TVh = 1-9 x 10 21 cm~ 2 x Ay (magnitude). 

The simulation of the CMD for the selected star- 
particles is shown in Fig. 2 where with a magnitude cut 
for true stars at above V > 20 mag about 8.8 x 10 4 star- 
particles within the galaxy coordinates defined above are 
retained. We recall that Fig. 2 shows absolute magnitude 
and luminosity, which are independent of the distance. To 
build the observational CMD we need to take into account 
distance and extinction for each star-particle, which are 
measured for each particle as described above. The resulting 
"observational" CMD is a shown in Fig. 5. With the IMF 
we assumed that the resulting CMD is representative of 
about 1.08 x 10 10 true stars. The difference in distance and 
extinction for each star-particle leads to a more smoothly 
distributed CMD. As a result, the groups of main sequence 
stars and red giant stars are hardly recognizable in Fig. 5. 
This is often observed in the photometric data of Galactic 
stars. 

7. Concluding remarks 

We first extended the theory of the stellar populations to 
include the phase-space description. We described the con- 
cepts of composite stellar population, simple stellar popu- 
lation, star formation rate, initial mass function, etc. us- 
ing the language of Statistical Mechanics. Now it is fair 




150 100 50 -50 -100 -150 



Fig. 4. Gas distribution in a Galactic coordinate of the simu- 
lated galaxy shown in Fig. 3. The region enclosed by the white 
lines are the selected region for constructing the CMD shown 
in Fig. 5. The colour coding is by density, so that bluer and 
more yellow points represent higher and lower density particles 
respectively. 




V-I 

Fig. 5. Observational CMD in apparent magnitudes and 
colours. The colour-coding is by density, so that bluer and more 
yellow points represent higher and lower density particles, re- 
spectively. 

to ask what we have gained from this rather formal ap- 
proach. These techniques have proven to be extremely pow- 
erful in other branches of physics e.g., for the study of re- 
laxation phenomena, electrical conduction, irreversible pro- 
cesses (see e.g., Landau & Lifshitz 2000) or in general for 
the development of the thermodynamics of non-equilibrium 
system (see e.g., de Groot & Mazur 1961) while have been 
more limited in the treatment of long-range forces e.g., 
the statistical mechanics of gravitational systems (see, e.g., 
Lynden-Bell & Wood 1968; Katz 2003), therefore it is im- 
portant to understand to what extent they can be applied 
to the stellar populations. The mutual benefit between 
Statistical Mechanics and stellar populations theory con- 
sists of the temporal evolution of the existence space we 
introduced, the projection of which onto the CMD is gov- 
erned by the fuel consumption theorem that drives the rel- 
ative number of stars in different evolutionary stages, hence 
cells of the CMD. Moreover, studying the coupling between 
dynamics and theory of the stellar populations imprinted 
in the existence space, may reveal unexplored connections. 
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The second achievement of this paper is indeed the ap- 
plication of this theory to develop a method for handling the 
synthetic CMDs that one would generate from N-body sim- 
ulations, which are to be compared with the CMDs of real 
galaxies when they are resolvable into stars. The method 
based on the concept of star frequencies per elemental cell 
of the CMD can easily simulate observational CMDs that 
contain huge numbers of stars. This study extends and com- 
pletes other similar studies in literature that dealt with 
synthetic CMDs such as those presented by in Carraro et 
al. (2001); Pasetto et al. (2012). This technique aims to 
interface N-body simulations to photometric observations 
and is developed with the perspective of the ever improv- 
ing capabilities of the N-body simulations in the future. 
Finally, thanks to its agility, the CMD tessellation method 
can be suitably interfaced with galaxy models based on star 
counts (see for instance Vallenari et al. 2006) and other hy- 
brid techniques based on the different approaches already 
present in literature such as the Galaxia model (see e.g., 
Sharma et al. 2011) or the Besancon model (see e.g. Robin 
et al. 2003) and it represents a natural extension of previous 
work (e.g. Tantalo et al. 2010). 

The tessellation code for generating the CMD for N- 
body simulations is available upon request from the au- 
thors. 

Acknowledgements. SP acknowledges D. Crnojevic and J. Hunt for 
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